%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% "Cyclical skill-biased technological change"
% Almut Balleer and Thijs van Rens
% Main file for estimation of shocks in a structural VAR framework
% This version for publication: October 2013
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% CHECK THE FOLLOWING AGAIN
% m-files and mat-files that you need in order to run this code:
% data-input: macrodata.m, Galidata.m, cpsdata_sa.m, DiCeccio_price.mat
% variables in the VAR: data_sec_pub
% estimation and identification: bvarest.m, ident.m, ident_sign.m,
% ident_sign_trijoint.m, ident_sign_tri.m
% other operations: VAR_lagtrend.m, vec.m, vec_to_mat.m, vardecomp.m, plots_pub.m, addresp.m
clear;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set options

% (1) specification
% determines the variables in the VAR and identification of shocks 
% Gali identification: 
% 1: productivity, hours
% 2: productivity, hours, premium, relative hours
% 3: productivity, hours, premium, price
% Fisher identification:
% 5: price, productivity, premium, hours
% SBT identification:
% 6: premium, productivity, hours
% 9: premium, productivity, hours, relative hours
% More additional variables
% 25: premium, relative hours, productivity, hours
% 26: premium, relative hours, productivity, hours, wage low skilled
% 27: premium, relative hours, productivity, hours, relative supply
% 28: premium, relative hours, productivity, hours, price
% 29: premium, relative hours, productivity, hours, price, investment
% 30: premium, relative hours, productivity, hours, price, investment,
% consumption
% 31: premium, relative hours, productivity, hours, price, investment,
% consumption, interest rate
% Joint identification:
% 32: price, premium, relative hours, productivity, hours
spec = 32;

if spec ==4 || spec == 7 || spec == 8 || spec >=10 && spec <=24
    display ('not a valid option, choose other spec')
    break
end

% (2) data
% different series for the skill premium and skill supply
prem_choice = 1; % 1: baseline series, 2: 'naive' measure
% different productivity series
prod_choice = 0; % set 0: BLS_data
% different series for the relative price of investment from DiCeccio
price_choice = 0; % 0: equipment price, 1: investment price
% different series for hours worked
hours_choice = 3; % set 0: BLS_data, 1: CPS private, 2: CPS all, 3: hours per capita

% (3) sample choice
% 1: short (1979:2-2000:4)
% 2: long (1979:2-2006:2)
% 5: 1979:2-2004:1
% 6: Gali (1948:1-1994:4)
% 7: Canova (1955:1-2000:4)
sample = 2;

if sample == 3 || sample == 4
    display ('not a valid option, choose other sample')
    break
end

% (4) identification
HORIZON_impresp = 20;  % number of periods in the impulse response functions 
HORIZON_ident = 80;    % identification horizon for long-run restrictions

% (5) specification of hours
% a) hours in levels or differences
do_level = 0;
% b) structural breaks in hours
do_break = 0;
break_choice = 0;
if do_break == 1
    % 0: dummy break  1: option with smooth HP   2: low-pass filter  
    % 3: deterministic trend, here 31: first order polynomial, 32: second
    % order polynomial, 33: third order polynomial
    break_choice = 2;
end

% (6) parameters for reduced form estimation
prior = 1;                  % type of prior, 1: Litterman, 2: Flat, 3: Kadiyala and Karlsson
nlagsvar = 8;               % lags in VAR
ndraws = 100;%1000;              % number of draws
phi = [0.2 .5 10^5];        % RATS manual default coefficients for tightness on own first lag
                            % other lags and exogenous variables
decay = 3;                  % harmonic decay of own lags, 1:linear decay
c = 1;

% (7) additional statistics
do_vardecomp = 0;

% (8) sign restrictions for SBT identification
do_sign = 2;

% select options until here, then run
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load data

% productivity
macrodata;
Galidata;
if prod_choice == 0
    if sample == 6 || sample == 7
      prod = log(gali_data(:,4));
    else
      prod = log(macro_data(:,8));
    end
end

% cps data on skill bias
cpsdata_sa;
if prem_choice == 1
    wage_premium = ourdata_sa(:,3);
    rel_empl = ourdata_sa(:,6);
elseif prem_choice == 2
    wage_premium = ourdata_rev_sa(:,6);
    rel_empl = ourdata_rev_sa(:,7);
end
rel_supply = ourdata_sa(:,8);
wagelow = ourdata_sa(:,10);

% relative price of investment
load('DiCeccio_price.mat', 'p_e08', 'p_i08', 'r_inv08');
if price_choice == 0
    if sample == 6 || sample == 7
        price = log(p_e08(5:238));
    else
        price = log(p_e08(129:238));
    end
elseif price_choice == 1
    price = log(p_i08(129:238));
end

% hours worked
if hours_choice == 0
    if sample == 6 || sample == 7
        hours = log(gali_data(:,3));
    else
        hours = log(macro_data(:,6)); 
    end
elseif hours_choice == 1
    hours = ourdata_sa(:,11);
elseif hours_choice == 2
    hours = ourdata_sa(:,12);
elseif hours_choice == 3
    if sample == 6 || sample == 7
        hours = log(gali_data(:,5));
    else
        hours = log(macro_data(:,11));
    end
end

% other macrodata
irate = macro_data(:,9);
consumption = log(macro_data(:,10));
investment = log(r_inv08(129:238));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% specification

% allowing for structural breaks
breaks = [];
if do_break == 1
    if break_choice == 0
        if sample == 1 || sample == 2 || sample == 5
            breaks = [72]; % 1997:I in differences
        elseif sample == 7
            breaks = [73 168];
        end
    end
end

% sample
rw = 0;

% define data and samplesize according to choices above
[data,nshocks,pndx,hndx,spndx,pricendx,wLndx,relendx,T,num_vars,frac_scalar] = data_spec_pub(...
    price,wage_premium,prod,hours,rel_supply,rel_empl,wagelow,irate,investment,consumption,spec,nlagsvar,prod_choice,do_break,break_choice,sample,do_level,breaks);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% estimating the BVAR

disp('Step I: Estimating the BVAR')
[Bet,Sigma,B,S] = bvarest(data,nlagsvar,hndx,c,ndraws,prior,T,phi,decay,do_level,breaks,break_choice);
% Bet are draws from posterior, B is the median of the posterior
% S is Sols, Sigma is Sols for each of the ndraws

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% structural identification

disp('Step II: Structural Identification')
if do_sign == 1 && spec == 9
    [Response,Respres]=ident_sign(Bet,Sigma,nlagsvar,ndraws,HORIZON_ident,c,...
        num_vars,nshocks,hndx,do_level);
    ndraws = size(Response,3);
elseif do_sign == 2
    if spec == 32
        [Response,Respres]=ident_sign_trijoint(Bet,Sigma,nlagsvar,ndraws,HORIZON_ident,c,...
            num_vars,nshocks,hndx,do_level);
        ndraws = size(Response,3);
    elseif spec >= 25
        [Response,Respres]=ident_sign_tri(Bet,Sigma,nlagsvar,ndraws,HORIZON_ident,c,...
            num_vars,nshocks,hndx,do_level);
        ndraws = size(Response,3);
    end
else
    [Response,Atilde,Respres]=ident(Bet,Sigma,nlagsvar,ndraws,HORIZON_ident,c,...
        num_vars,nshocks,hndx,do_level);
end

% additional responses for output, wages of low and high skilled and
% employment for low and high skilled
[num_vars_ext,Respext,Respextres,yndx,wHndx,hHndx,hLndx] = ...
    addresp(Response,Respres,spec,num_vars,ndraws,pndx,hndx,spndx,wLndx,relendx,frac_scalar);

% CHECK THIS AGAIN
if spec == 5
        save(['Response_spec' int2str(spec) '_sample' int2str(sample) '_prem_' int2str(prem_choice) '_prod' int2str(prod_choice) '_price' int2str(price_choice) '_hours' int2str(hours_choice) '_break' int2str(do_break) int2str(break_choice) '_est' int2str(prior) int2str(nlagsvar) int2str(decay)],'Response')
end
% sorting the Responses
Respsort = zeros(size(Respext));
for j=1:(size(Respext,1))
    Respsort(j,:,:) = sort(Respext(j,:,:),3);
end;
Respressort = zeros(size(Respextres));
for j=1:size(Respres,1)
    Respressort(j,:,:) = sort(Respextres(j,:,:),3);
end;

% Median and 68% credible set
Respmed(:,:) = Respsort(:,1:HORIZON_impresp,fix(0.5*ndraws));
Resplow(:,:) = Respsort(:,1:HORIZON_impresp,fix(0.16*ndraws));
Resphigh(:,:) = Respsort(:,1:HORIZON_impresp,fix(0.84*ndraws));
Respresmed(:,:) = Respressort(:,1:HORIZON_impresp,fix(0.5*ndraws));
Respreslow(:,:) = Respressort(:,1:HORIZON_impresp,fix(0.16*ndraws));
Respreshigh(:,:) = Respressort(:,1:HORIZON_impresp,fix(0.84*ndraws));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% variance decomposition

if do_vardecomp == 1
    disp('Step III: Variance Decomposition')
    % calculate the variance decomposition from the impulse responses
    [frac_rev_med,frac_rev_low,frac_rev_high] = vardecomp(Respext,num_vars,num_vars_ext,HORIZON_ident,ndraws);
    
    % display for output, hours, premium, productivity
    output_mat = [];
    hours_mat = [];
    premium_mat = [];
    prod_mat = [];
    rele_mat = [];
    for n=1:num_vars
        output_mat = [output_mat
            100*frac_rev_med(yndx,n,1),0,100*frac_rev_med(yndx,n,8),0,100*frac_rev_med(yndx,n,16),0,100*frac_rev_med(yndx,n,32),0;
            100*frac_rev_low(yndx,n,1),100*frac_rev_high(yndx,n,1),100*frac_rev_low(yndx,n,8),100*frac_rev_high(yndx,n,8),100*frac_rev_low(yndx,n,16),100*frac_rev_high(yndx,n,16),100*frac_rev_low(yndx,n,32),100*frac_rev_high(yndx,n,32)];
        hours_mat = [hours_mat
            100*frac_rev_med(hndx,n,1),0,100*frac_rev_med(hndx,n,8),0,100*frac_rev_med(hndx,n,16),0,100*frac_rev_med(hndx,n,32),0;
            100*frac_rev_low(hndx,n,1),100*frac_rev_high(hndx,n,1),100*frac_rev_low(hndx,n,8),100*frac_rev_high(hndx,n,8),100*frac_rev_low(hndx,n,16),100*frac_rev_high(hndx,n,16),100*frac_rev_low(hndx,n,32),100*frac_rev_high(hndx,n,32)];
        premium_mat = [premium_mat
            100*frac_rev_med(spndx,n,1),0,100*frac_rev_med(spndx,n,8),0,100*frac_rev_med(spndx,n,16),0,100*frac_rev_med(spndx,n,32),0;
            100*frac_rev_low(spndx,n,1),100*frac_rev_high(spndx,n,1),100*frac_rev_low(spndx,n,8),100*frac_rev_high(spndx,n,8),100*frac_rev_low(spndx,n,16),100*frac_rev_high(spndx,n,16),100*frac_rev_low(spndx,n,32),100*frac_rev_high(spndx,n,32)];
        prod_mat = [prod_mat
            100*frac_rev_med(pndx,n,1),0,100*frac_rev_med(pndx,n,8),0,100*frac_rev_med(pndx,n,16),0,100*frac_rev_med(pndx,n,32),0;
            100*frac_rev_low(pndx,n,1),100*frac_rev_high(pndx,n,1),100*frac_rev_low(pndx,n,8),100*frac_rev_high(pndx,n,8),100*frac_rev_low(pndx,n,16),100*frac_rev_high(pndx,n,16),100*frac_rev_low(pndx,n,32),100*frac_rev_high(pndx,n,32)];
        rele_mat = [rele_mat
            100*frac_rev_med(relendx,n,1),0,100*frac_rev_med(relendx,n,8),0,100*frac_rev_med(relendx,n,16),0,100*frac_rev_med(relendx,n,32),0;
            100*frac_rev_low(relendx,n,1),100*frac_rev_high(relendx,n,1),100*frac_rev_low(relendx,n,8),100*frac_rev_high(relendx,n,8),100*frac_rev_low(relendx,n,16),100*frac_rev_high(relendx,n,16),100*frac_rev_low(relendx,n,32),100*frac_rev_high(relendx,n,32)];
    end
    display(output_mat);
    display(hours_mat);
    display(premium_mat);
    display(prod_mat);
    display(rele_mat);
else
    disp('Step III: Variance Decomposition not performed')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plots

plots_pub(Respmed,Resplow,Resphigh,spec,HORIZON_impresp,num_vars_ext,sample,rw,prem_choice,prod_choice,price_choice,hours_choice,do_break,break_choice,...
    prior,nlagsvar,decay,pndx,hndx,spndx,pricendx,wLndx,relendx,wHndx,hHndx,hLndx,do_level,do_sign)
